package edu.amm.fanal.methods;

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

import org.apache.log4j.Logger;

import edu.amm.algorithms.Cholesky;
import edu.amm.api.Bounds;
import edu.amm.fanal.api.BoundaryProblem;
import edu.amm.fanal.api.BoundsTransformation;
import edu.amm.fanal.api.Coordinates;
import edu.amm.fanal.api.Coordinates.CoordinatesType;
import edu.amm.fanal.api.Function;
import edu.amm.fanal.api.FunctionalSpace;
import edu.amm.fanal.api.FunctionalSpace.SpaceType;
import edu.amm.fanal.api.LinearProblem;
import edu.amm.fanal.api.Transformation;
import edu.amm.fanal.api.ZeroTransformation;

/**
 * Метод моментов для решения краевой задачи для
 * линейного дифференцирального уравнения второго порядка.
 * 
 * @author Иванов Илья
 * @since 2013-05-10
 *
 */
public abstract class MomentsMethod extends Method<LinearProblem> {
	
	/** Задача с нулевыми краевыми условиями, которую
	 * необходимо решить конкретному методу моментов */
	protected LinearProblem zeroProblem;
	
	/** Пространство, в котором работает конкретный метод моментов. */
	protected FunctionalSpace space;
	
	/** Базис для построения текущего решения. */
	protected Function[] basis;
	
	private SpaceType spaceType;
	private double epsilon = DEFAULT_ACCURACY;
	
	private CoordinatesType coordinatesType;
	
	/**
	 * @param spaceType Тип пространства, в котором работает метод.
	 * @param coordinatesType Тип координатной системы.
	 */
	public MomentsMethod(SpaceType spaceType, CoordinatesType coordinatesType) {
		if (spaceType == null)
			throw new RuntimeException("Не задано пространство, в котором работает метод моментов");
		if (coordinatesType == null)
			throw new RuntimeException("Не задан тип координатной системы для метода моментов");
		
		this.spaceType = spaceType;
		this.coordinatesType = coordinatesType;
	}
	
	/**
	 * @param epsilon Требуемая норма разности двух последних
	 * приближений к решению задачи. Значение по умолчанию
	 * равно {@link #DEFAULT_ACCURACY}
	 */
	public void setAccuracy(double epsilon) {
		if (epsilon < 0)
			throw new RuntimeException("Близость по норме должна быть неотрицательной");
		
		this.epsilon = epsilon;
	}
	
	/**
	 * @return Название реализации метода моментов.
	 */
	public abstract String getName();
	
	/**
	 * Коэффициент матрицы СЛАУ для получения текущего приближения к
	 * решению задачи.
	 * 
	 * @param i Номер строки матрицы, начиная с 1.
	 * @param j Номер столбца матрицы, начиная с 1.
	 * @return Коэффициент (i, j).
	 */
	protected abstract double getMatrixCoefficient(int i, int j);
	
	/**
	 * Элемент вектора правой части СЛАУ для получения текущего приближения к
	 * решению задачи.
	 * 
	 * @param i Номер элемента, начиная с 1.
	 * @return <code>i</code>-й элемент вектора правой части СЛАУ.
	 */
	protected abstract double getRightHandElement(int i);
	
	protected Function solveProblem(LinearProblem problem) {
		log.info("Метод моментов. Реализация: " + getName());
		
		List<Transformation> transformations = new ArrayList<>();
		
		Transformation transformation = new ZeroTransformation(problem);
		LinearProblem targetProblem = transformation.getTransformed();
		
		log.info("Применено преобразование линейной задачи к задаче с нулевыми краевыми условиями");
		
		transformations.add(transformation);
		
		if (coordinatesType == CoordinatesType.INTEGRAL_LEGENDRE) {
			Bounds wasBounds = targetProblem;
			transformation = new BoundsTransformation(targetProblem, Bounds.LEGENDRE_BOUNDS);
			targetProblem = transformation.getTransformed();
			
			log.info(String.format("Применено преобразование границ линейной задачи с отрезка [%f, %f] к отрезку [%f, %f]",
				wasBounds.getLeftBound(), wasBounds.getRightBound(), targetProblem.getLeftBound(), targetProblem.getRightBound()));
			
			transformations.add(transformation);
		}
		
		zeroProblem = targetProblem;
		Function solution = solveZeroProblem();
		
		Collections.reverse(transformations);
		for (Transformation transform : transformations)
			solution = transform.restoreSolution(solution);
		
		return solution;
	}
	
	private Function solveZeroProblem() {
		log.info("Решение задачи с нулевыми краевыми условиями");
		
		space = spaceType.instantiate(zeroProblem);
		
		log.info("Используется функицональное пространство: " + space.getSpaceName());
		
		Coordinates coordinates = coordinatesType.instantiate(zeroProblem);
		
		log.info("Используется координатная система: " + coordinates.getName());
		
		Function solution = null;
		double accuracy = Double.MAX_VALUE;
		int subspaceDim = FIRST_SUBSPACE_DIM;
		
		while (accuracy > epsilon && subspaceDim <= MAX_SUBSPACE_DIM) {
			log.info(String.format("Итерация: %d. Размерность подпространства: %d",
				subspaceDim - FIRST_SUBSPACE_DIM + 1, subspaceDim));
			
			basis = coordinates.getCoordinates(subspaceDim);
			
			Function previous = solution;
			
			log.info("Вычисление матрицы СЛАУ...");
			
			double[][] matrix = getMatrix(subspaceDim);
			
			log.info("Вычисление правой части СЛАУ...");
			
			double[] rightHand = getRightHand(subspaceDim);
			
			log.info("Решение СЛАУ методом Холецкого...");
			
			double[] coefficients = new Cholesky().solve(matrix, rightHand);
			
			solution = buildZeroSolution(coefficients);
			
			log.info(String.format("Получено приближённое решение задачи в подпространстве размерности %d", subspaceDim));
			
			if (previous != null) {
				accuracy = getAccuracy(solution, previous);
				
				log.info(String.format("Разность по норме с решением, полученным на предыдущей итерации: %.17f", accuracy));
			}
			
			subspaceDim++;
		}
		
		if (accuracy > epsilon)
			log.warn("Процесс решения методом моментов остановлен по превышении максимально допустимой размерности подпространства");
		else
			log.info("Требуемая точность в методе моментов достигнута");
		
		return solution;
	}
	
	private double[][] getMatrix(int subspaceDim) {
		double[][] matrix = new double[subspaceDim][subspaceDim];
		for (int i = 0; i < matrix.length; i++)
			for (int j = 0; j < matrix[i].length; j++)
				matrix[i][j] = getMatrixCoefficient(i, j);
		
		return matrix;
	}
	
	private double[] getRightHand(int subspaceDim) {
		double[] rightHand = new double[subspaceDim];
		for (int i = 0; i < rightHand.length; i++)
			rightHand[i] = getRightHandElement(i);
		
		return rightHand;
	}
	
	private Function buildZeroSolution(final double[] coefficients) {
		final Function[] basis = this.basis;
		return new Function(BoundaryProblem.SOLUTION_DIM) {
			
			private Function derivative = new Function(getDim()) {
				
				private Function[] dbasis = new Function[basis.length];
				
				{
					for (int i = 0; i < dbasis.length; i++)
						dbasis[i] = basis[i].differentiate(0);
				}
				
				protected double getValue(double... x) {
					double value = 0;
					for (int i = 0; i < coefficients.length; i++)
						value += coefficients[i] * dbasis[i].evaluate(x);
					
					return value;
				}
			};
			
			protected double getValue(double... x) {
				double value = 0;
				for (int i = 0; i < coefficients.length; i++)
					value += coefficients[i] * basis[i].evaluate(x);
				
				if (value > THRESHOLD)
					log.error(value);
				return value;
			}
			
			protected Function getDerivative(int variableNumber) {
				return derivative;
			}
		};
	}
	
	private double getAccuracy(final Function solution, final Function previous) {
		return space.norm(Function.subtraction(solution, previous));
	}
	
	/** По умолчанию близость по норме двух последних
	 * приближений к решению задачи. */
	public static final double DEFAULT_ACCURACY = 1e-3;
	
	/** Размерность подпространства, в котором лежит первое
	 * приближение к решению задачи. */
	public static final int FIRST_SUBSPACE_DIM = 5;
	
	/** Максимальная размерность подпространтсва, в котором
	 * может лежать приближение к решению задачи */
	private static final int MAX_SUBSPACE_DIM = 50;
	
	private static final double THRESHOLD = 1e+9;
	private static Logger log = Logger.getLogger(MomentsMethod.class);
}